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Abstract 

Given a zero-dimensional ideal / c K[jfi ,... ,x„] of degree D, the trans- 
formation of the ordering of its Grobner basis from DRL to LEX is a key 
f) , step in polynomial system solving and turns out to be the bottleneck of the 

ry*\ • whole solving process. Thus it is of crucial importance to design efficient 

• , algorithms to perform the change of ordering. 

Q ■ The main contributions of this paper are several efficient methods for the 

change of ordering which take advantage of the sparsity of multiplication 
matrices in the classical FGLM algorithm. Combing all these methods, we 
propose a deterministic top-level algorithm that automatically detects which 
method to use depending on the input. As a by-product, we have a fast im- 
C") \ plementation that is able to handle ideals of degree over 40000. Such an 

implementation outperforms the Magma and Singular ones, as shown by our 
experiments. 

First for the shape position case, two methods are designed based on the 
Wiedemann algorithm: the first is probabilistic and its complexity to com- 
plete the change of ordering is 0{D(N\ + nlog(D))), where N\ is the num- 
ber of nonzero entries of a multiplication matrix; the other is deterministic 
and computes the LEX Grobner basis of \fl via Chinese Remainder The- 
orem. Then for the general case, the designed method is characterized by 
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H ■ the Berlekamp-Massey-Sakata algorithm from Coding Theory to handle the 

multi-dimensional linearly recurring relations. Complexity analyses of all 
proposed methods are also provided. 

Furthermore, for generic polynomial systems, we present an explicit for- 
mula for the estimation of the sparsity of one main multiplication matrix, 
and prove its construction is free. With the asymptotic analysis of such spar- 
sity, we are able to show for generic systems the complexity above becomes 
0{^6jmtD 2 + ,J ir). 
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Email addresses: Jean-Charles.Faugere@inria.fr (Jean-Charles Faugere), Chenqi.Mou@lip6.fr 
(Chenqi Mou) 



1 Introduction 

1.1 Motivation 

Grobner basis is an important tool in computational ideal theory ||9] [T3j [6), es- 
pecially for polynomial system solving. For a given ideal and term ordering, the 
Grobner basis of this ideal with respect to (w.r.t.) the term ordering is a set of gen- 
erators with good properties, such that manipulation of the ideal can be achieved 
with these generators. 

The term ordering plays an important role in the theory of Grobner bases. It is 
well-known that Grobner bases w.r.t. different term orderings are also different and 
possess different theoretical and computational properties. For example, Grobner 
bases w.r.t. the lexicographical ordering (LEX) have good algebraic structures and 
are convenient to use for polynomial system solving, while those w.r.t. the degree 
reverse lexicographical ordering (DRL) are computationally easy to obtain. There- 
fore, the common strategy to solve a polynomial system is to first compute the 
Grobner basis of the ideal defined by the system w.r.t. DRL, change its ordering to 
LEX, and perhaps further convert the LEX Grobner basis to triangular sets [24] or 
Rational Univariate Representation [30|. That is one of the main usages of algo- 
rithms for the change of ordering. 

However, the computation of Grobner bases greatly enhanced recently |[T7l[T8l . 
the step to change the ordering of Grobner bases has become the bottleneck of 
the whole solving process (see Section |7J. Hence it is of crucial significance to 
design efficient algorithms for the change of ordering. The purpose of this paper is 
precisely to provide such efficient algorithms. 

Furthermore, some practical problems can be directly modeled as the change 
of ordering of Grobner bases. For example, the Grobner basis of an ideal derived 
from the AES-128 cryptosystem w.r.t. a certain term ordering (other than LEX) has 
been obtained [ 10], and it may lead to a successful cryptanalysis on this system if 
one is able to convert its term ordering to LEX. And the decoding of some cyclic 
codes can also be regarded as a problem of changing the term ordering ll25l . 

1.2 Related works 

Several algorithms for the change of ordering have already existed, for example 
the FGLM algorithm for the zero-dimensional case [19] and the Grobner walk for 
the generic case |[T2l . Similar algorithms have also been proposed to change the 
orderings of triangular sets |]29] H4l or using the LLL algorithm [4 1 in the bivariate 
case. 

Among them, the FGLM algorithm, only applicable to the zero-dimensional 
case, is an efficient one. The number of field operations it needs to complete the 
change of ordering is 0(nD 3 ), where n is the number of variables, and D is the 
degree of the given ideal / C K[xi, . . . ,x n ]. Its efficiency may be due to the fact that 
it reduces the problem of change of ordering to linear algebra operations. Such 



a connection is achieved through the multiplication matrix 7} (i = l,...,n) used 
in this algorithm, which represents the multiplication by x, in the quotient ring 
K[jci ,... ,x n ]/I viewed as a vector space. These matrices are sparse, even when the 
input polynomial system is dense (see Section©. And in this paper we take advan- 
tage of this sparsity structure to obtain fast FGLM algorithms with good complexity 
and performances. 

1.3 Our contributions 

We first study the particular but important case when the zero -dimensional ideal I is 
in shape position. Two methods based on the Wiedemann algorithm are proposed 
to compute the Grobner bases of / or \fl w.r.t. LEX. They both make use of the 
sparsity by constructing the linearly recurring sequence 

[<r,r/e):/ = 0,...,2D-l], 

where r is a vector and e = ( 1 , 0, . . . )' is the vector representing 1 in K[x\ , . . . ,x n ] /I. 
It is easy to see that the minimal polynomial f\ in K[jei] of this linearly recurring 
sequence is indeed a polynomial in the Grobner basis of I w.r.t LEX (x\ < ■■■ <x n ) 
when deg(/i) = D, and it can be computed by applying the Berlekamp-Massey 
algorithm [36]. Furthermore, we show how to recover efficiently the other polyno- 
mials in the Grobner basis by solving structured (Hankel) linear systems. Hence, 
we are able to complete the first method for the change of ordering to LEX for 
ideals in shape position with complexity 0{D(N\ +wlog(D))), where Ni is the 
number of nonzero entries in T\ . When n <^D this almost matches the complexity 
of computing the minimal polynomial. 

The other method for the shape position case uses the deterministic Wiede- 
mann algorithm, which can always return the correct univariate polynomial in the 
Grobner basis w.r.t. LEX. Making use of the Chinese Remainder Theorem, this 
method adapts and extends the previous one to recover the Grobner basis of y/l, 
instead of I. Thus it is suitable to those problems where the zeros, instead of 
the multiplicities, are of interest. For ideals in shape position, this deterministic 
method can always return the Grobner basis of \fl w.r.t. LEX with the complexity 
0(D(Ni +D(log(D)loglog(D) + «))). 

We also briefly discuss how to apply an incremental variant of the Wiedemann 
algorithm to compute the univariate polynomial, which is of special importance 
among all the polynomials in the Grobner basis. Such an variant has a complexity 
sensitive to the output, namely the degree of the univariate polynomial, and is 
efficient when this degree is small. 

Then for general ideals to which the methods above may be no longer applica- 
ble, we follow the idea above by generalizing the linearly recurring sequence to a 
^-dimensional mapping 

E:(s u ...,s n )^^(r,T( 1 ---T^e). 



The minimal set of generating polynomials (w.r.t. a term ordering) for the linearly 
recurring relation determined by E is essentially the Grobner basis of the ideal 
defined by E, and this polynomial set can be obtained via the Berlekamp-Massey- 
Sakata (BMS for short hereafter) algorithm from Coding Theory [f32l 1331 . With 
modifications of the BMS algorithm, we design a method to change the ordering 
in the general case. The complexity of this algorithm is 0(nD(N + NND)), where 
N is the maximal number of nonzero entries in matrices T\ , . . . , T„, while N and N 
are respectively the number of polynomials and the maximal term number of all 
polynomials in the resulting Grobner basis. 

Combing all these methods above, we present a deterministic top-level algo- 
rithm, which is able to choose automatically which method to use according to 
the input. The efficiency of the proposed methods is verified by experiments. The 
current implementation outperforms those of FGLM in Magma and Singular. Take 
a randomly generated quadratic polynomial system of 13 variables for example, it 
generates an ideal in shape position of degree 8192. For such an ideal, the change 
of ordering to LEX can be achieved in 193.5 seconds: this is 54 times faster than 
the corresponding Magma function. As shown in Table|7J zero-dimensional ideals 
over a prime field of degree greater than 40000 are now tractable. 

Furthermore, the performances of these methods are heavily dependent on the 
sparsity of the multiplication matrices, especially T\ for the shape position case. 
In general we assume the multiplication matrices known. However, for generic 
polynomial systems consisting of n-variate polynomials of degree d, the sparsity 
of T\ is investigated, and we are able to give an explicit formula to compute the 
number of dense columns in T\ and show indeed its construction is free. These 
results furnish a complete complexity analysis of the proposed method for generic 
polynomial systems. Then with an asymptotic analysis of the number of dense 
columns as d tends to +oo, we show the complexity of the first method for the shape 
position case becomes 0(y // 6/nnD 2+ ~^~) for generic systems. Such simplified 
complexity is better than that of FGLM with smaller constant and exponent. 

1.4 What is new 

To be self-contained, this paper also includes results obtained in ifTSI in a refined 
way for description. However, several original extensions have also been presented 
here, making the discussion on this subject more comprehensive: (1) For ideals in 
shape position, one new algorithm is proposed based on the deterministic Wiede- 
mann algorithm. Compared with the previous probabilistic one, this algorithm 
becomes deterministic and aims at the Grobner basis of the radical of the input 
ideal. (2) The multiplication matrices are assumed known in [15], but here for the 
multiplication matrix 7\ which is of special importance, its sparsity, together with 
the asymptotic behaviors, and construction cost are analyzed for generic polyno- 
mial systems. Such a study furnishes a complete understanding for the complexity 
of the change of ordering for generic systems, with construction of multiplication 
matrices also considered. (3) The proof of Theorem 14- 1 1 is further simplified via 



introduction of known results in the literature. 

1.5 Paper structure 

The organization of this paper is as follows. Related preparatory algorithms used in 
this paper, along with some notations, are first reviewed in Section|2] Then Section 
[3] is devoted to the shape position case, where two methods with their complexity 
analyses are exploited. The method based on the BMS algorithm for the general 
case is presented in Section |4] Section [5] combines all the previous methods to a 
top-level algorithm. The sparsity of T\ is studied in Section [6] and experimental 
results are provided in Section [7] 

2 Backgrounds: FGLM and BMS algorithms 

Let K[x\, . . . ,x n ] be the rc-variate polynomial ring over a field K, with variables 
ordered as x\ < • • • < x n . Suppose G\ is the Grobner basis of a O-dimensional ideal 
I C K[xi, . . . ,x n ] w.r.t. a term ordering <i. Given another term ordering <2, one 
wants to compute the Grobner basis G2 of I w.r.t. it. Denote by D the degree of 
I, that is, the dimension of K[xi,. . . ,x n ]/I as a vector space. These notations are 
fixed hereafter in this paper. 

2.1 FGLM algorithm 

The FGLM algorithm is one to perform the change of ordering of Grobner bases of 
O-dimensional ideals efficiently [19]. The reason why it is fast may be due to the 
idea that it reduces the problem of ordering change to linear algebra operations in 
the quotient ring K[xi, . . . ,x n ]/I. Such a reduction is realized in the following way. 

First one computes the canonical basis of K[xi,...,x n ]/(Gi) and orders its 
elements according to <i. Let B = [ei,...,Ed] be the ordered basis. Then £1 
will always equal 1, for <i is a term ordering. Given a variable x,-, for each el- 
ement Sj in B, one can compute the normal form of EjXi w.r.t. G\, denoted by 
NormalForm(£ 7 o; r ). This normal form, viewed as an element of K[xi , . . . ,x n ]/(G\), 
can be further written as a linear combination of B. Writing the coefficients as a 
column vector, one can construct a D x D matrix 7} by adjoining all the column 
vectors for j = 1, . . . ,D. This matrix is called the multiplication matrix of x\. It is 
not hard to verify that all 7} commute: 7}7) = 7)7} for i,j = 1, . . . ,n. 

Next one handles all the terms in K[jci,. . . ,x n ] one by one following <2. For 
each term x s with s = (si,... ,s n ), its coordinate vector w.r.t. B can be computed 
by 

where e = (1,0, .. . ,0) f is the coordinate vector of 1. Then criteria proposed in 
FGLM guarantee that once a linear dependency of the coordinate vectors of com- 



puted terms 

£c s u s =0 (1) 

seS 

is found, a polynomial /6G2 can be directly derived in the following form 

f = x l + £ ^ s , (2) 

seS^l Cl 

where x l is the leading term of / w.r.t. <2 (denoted by lt(/)) [ 19]. 

As can be seen now, all one needs to do to obtain the Grobner basis G2 is 
computing the coordinate vector of each term one by one, and checking whether a 
linear dependency of these vectors occurs after a new vector is computed, which 
can be realized by maintaining an echelon form of the matrix whose columns are 
coordinate vectors of previously computed terms. These steps are merely matrix 
manipulations from linear algebra. A trivial upper bound for the number of terms 
to consider is D + 1 because of the vector size. 

2.2 BMS algorithm 

The BMS algorithm from Coding Theory is a decoding algorithm to find the gen- 
erating set of the error locator ideal in algebraic geometry codes |[32l]33l]3TI . From 
a more mathematical point of view, it computes the set of minimal polynomials 
(w.r.t. a term ordering <) of a linearly recurring relation generated by a given 
multi-dimensional array. It is a generalization of the Berlekamp-Massey algo- 
rithm, which is applied to Reed-Solomon codes to find the generating error locator 
polynomial, or mathematically the minimal polynomial of a linearly recurring se- 
quence. 

The BMS algorithm, without much modification, can also be extended to a 
more general setting of order domains lTT3l I2T1 . Combining with the Feng-Rao 
majority voting algorithm BOl . this algorithm can often decode codes with more 
with ( dmm — 1)/2 errors if the error locations are general [T , where d m [ n is the min- 
imal distance. Next a concise description of the BMS algorithm is given, focusing 
on its mathematical meanings. 

As a vector u = (u\,...,u n ) G Z> and a term x u = x" 1 ■■■x" n E K[xi,.. . ,x n ] 
are 1-1 corresponding, usually we do not distinguish one from the other. A map- 
ping E : Z> — ► K is called a n-dimensional array. In Coding Theory, the array 
E is usually a syndrome array determined by the error word [31 1 . Besides the term 
ordering, we define the following partial ordering: for two terms u = (u\,...,u n ) 
and v = (vi, .. . ,v n ), we say that u -< v if iij < v, for i = 1, .. . ,n. 

Definition 2.1 Given a polynomial / = £ s / s ^ s £ &[xi,. • • ,x n ], a ^-dimensional 
mapping E is said to satisfy the n-dimensional linearly recurring relation with 
characteristic polynomial f if 

J £fsE s+r = 0, Vr^O. (3) 



The set of all characteristic polynomials of the ^-dimensional linearly recurring 
relation for the array E forms an ideal, denoted by 1(E). Again in the setting of 
decoding when E is a syndrome array, this ideal is called the error locator ideal 
for E, and its elements are called error locators. The definition of 1(E) used here 
in this paper follows 11311 . and one can easily see that this definition is equivalent 
to that in US by EB Thereom 23]. 

Furthermore, the set of minimal polynomials for 1(E) w.r.t. <, which the BMS 
algorithm computes, is actually the Grobner basis of 1(E) w.r.t. < 031 Lemma 5], 
The canonical basis of K[x\ ,... ,x n ]/I(E) is also called the delta set of E, denoted 
by A(is). The term "delta set" comes from the property that if u G Z> is contained 
in A (is), then A (is) also contains all elements v G Z> such that v -< u. 

Instead of studying the infinite array is as a whole, the BMS algorithm deals 

with a truncated subarray of E up to some term u according to the given term 

ordering <. A polynomial / with lt(/) = s is said to be valid for E up to u if 

either u )f s or 

Y,ftE t +r = 0, W(0-<r<u-s). 
t 

E may be omitted if no ambiguity occurs. A polynomial set is said to be valid up 
to u if each its polynomial is so. 

Similarly to FGLM, the BMS algorithm also handles terms in K[xi,. . . ,x n ] one 
by one according to <, so that the polynomial set F it maintains is valid up to the 
new term. Suppose F is valid up to some term u. When the next term of u w.r.t. 
<, denoted by Next (it), is considered, the BMS algorithm will update F so that 
it keeps valid up to Next(u). Meanwhile, terms determined by Next(-u) are also 
tested whether they are members of A (is). Therefore, more and more terms will be 
verified in A(is ) as the BMS algorithm proceeds. The set of verified terms in A(is ) 
after the term u is called the delta set up to u and denoted by A(u). Then we have 

A(l) C • • • C A(u) C A(Next(u)) C • • • C A(E). 

After a certain number of terms are considered, F and A (it) will grow to the 
Grobner basis of 1(E) and A(E) respectively. 

Next only the outlines of the update procedure mentioned above, which is also 
the main part of the BMS algorithm, are presented as Algorithm [TJfor convenience 
of later use. More details will also be provided in Section |4] One may refer to 
lOTl [TBI for a detailed description. In Algorithm [TJ below, the polynomial set G, 
called the witness set, is auxiliary and will not be returned with F in the end of the 
BMS algorithm. 

3 Shape position case: probabilistic, deterministic and in- 
cremental algorithms 

In this section, the case when the ideal I is in shape position is studied. 



Algorithm 1: (F + ,G + ) := BMSUpdate(F,G,Next(w),£) 
Input: 

F, a minimal polynomial set valid up to it; 

G, a witness set up to u; 
Next(tt), a term; 

E, a ^-dimensional array up to Next (it). 
Output: 

F + , a minimal polynomial set valid up to Next(it); 
G + , a witness set up to Next(it). 

1. Test whether every polynomial in F is valid up to Next (it) 

2. Update G + and compute the new delta set up to Next(it) accordingly 

3. Construct new polynomials in F + such that they are valid up to Next(it) 



Definition 3.1 An ideal / C IK[jci,...,jc„] is said to be in shape position if its 
Grobner basis w.r.t LEX is of the following form 

[fi{xi),X2-f 2 (xi),...,x n -f n (xi)]. (4) 

One may easily see that / here is O-dimensional and deg(/i) = D. 

Such ideals take a large proportion in all the consistent ideals and have been 
well studied and applied [5, 30]. The special structure of their Grobner bases en- 
ables us to design specific and efficient methods to change the term ordering to 
LEX. In the following, methods designed for different puiposes, along with their 
complexity analyses, are exploited. 

Throughout this section, we assume the multiplication matrix T\ is nonsingular. 
Otherwise, one knows by the Stichelberger's theorem (cf. [30, Theorem 2.1]) that 
x\ = will be a root of the univariate polynomial in /'s Grobner basis w.r.t. LEX, 
and sometimes the polynomial system can be further simplified. 

3.1 Probabilistic algorithm to compute Grobner basis of the ideal 

3.1.1 Algorithm description 

Given a O-dimensional ideal /, if the univariate polynomial f\ (x\ ) in its Grobner 
basis w.r.t. LEX is of degree D, then we know / is in shape position. 

The way to compute such a univariate polynomial is the Wiedemann algorithm. 
Consider now the following linearly recurring sequence 

s=[(r,T{e):i = 0,...,2D-l], (5) 

where r is a randomly generated vector in K( Dxl \ T\ is the multiplication matrix 
of x\ , e is the coordinate vector of 1 w.r.t the canonical basis of K[xi , . . . ,x„]/I, and 



(•, •) takes the inner product of two vectors. It is not hard to see that the minimal 
polynomial f\ of the sequence s is a factor of f\. As D is always a bound on the 
size of the linearly recurring sequence, the Berlekamp-Massey algorithm can be 
applied to the sequence s to compute f\ . Furthermore, if deg (f\ ) = D, then f\ = f\ 
and / can be verified in shape position. 

Suppose deg(/i) = D holds and /; in (@| is of the form /, = Y^=o c i,kx\ for 
i = 2, ... ,n. Then computing the whole Grobner basis of / w.r.t. LEX reduces 
to determining all the unknown coefficients c,-^. Before we show how to recover 
them, some basic results about linearly recurring sequences are recalled. 

Definition 3.2 Let s = [so,s\,S2,- ••] be a sequence of elements in IK and d an 
integer. The d x d Hankel matrix is defined as 



H d (s) 






S2 



S2 
S3 



Sd-1 s d Sd+\ 



Theorem 3.1 (|22]) Let s 



-d 



Sd-1 
Sd 

Sld-2 



be a linearly recurring sequence. Then 



the minimal polynomial M^ s '(x) = Y^i=Q m i x '°f tne sequence s is such that: 
(i) d = rank(Hd(s)) = rank (H((s)) for all i > d; 
(ii) ker(//^ + i (s)) is a vector space of dimension 1 generated by (mo, mi , . 



,m d ) 



For each i = 2, . . . , n, as xi — L^Lq 1 c i,k x \ £ I> one has NormalForm (x,- — L^Lq 1 c i,k x \ ] 



0, thus 



D-l 



T ie = Y, cik ■ rfe. 



k=0 



Multiplying T( and taking the inner product with a random vector r to both hands 
for j = 1, . . . ,D — 1, one can further construct D linear equations 

D-l 

(r, T(vi) = £ c iik ■ (r, T? +i e) , j = 0, . . . ,D - 1 . 



(6) 



k=0 



k +J. 



With aji considered as unknowns, the coefficient matrix H with entries (r, Tj e 



is indeed a D x D Hankel matrix, and thus invertible by Theorem 13. II Furthermore, 
the linear equation set © with the Hankel matrix H can be efficiently solved JH. 
All the solutions of these linear systems for i = 2, . . . ,n will lead to the Grobner 
basis we want to compute. 

The method above is summarized in the following algorithm, whose termina- 
tion and correctness are direct results based on previous discussions. The sub- 
function BerlekampMassey() is the Berlekamp-Massey algorithm, which takes a 
sequence over K as input and returns the minimal polynomial of this sequence [36]. 

Remark 3.1 As can be seen from the description of Algorithm |2] such a method 
is a probabilistic one. That is to say, it can return the correct Grobner basis w.r.t. 
LEX with probabilities, and may also fail even when / is indeed in shape position. 



Algorithm 2: Shape position (probabilistic) G 2 '■= ShapePro(Gi, <i) 
Input: Gi , Grobner basis of a O-dimensional ideal / C K.[x\ , . . . ,x n ] w.r.t. < 1 
Output: G 2 , Grobner basis of / w.r.t. LEX if the polynomial returned by 
BerlekampMassey() is of degree D; Fail, otherwise. 

Compute the canonical basis of K[xi, . . . ,x n ]/(Gi) and multiplication 
matrices T\,...,T n ; 
e:=(l,0,...,0)' £l( Dx1 '; 
Choose vq = r G IK> Dxl ) randomly; 
forz = 1,...,2D-1 do 

I n := (r/)n_ i; 

end 

Generate the sequence s := [(r,, e) : 2 = 0, . . . , 2D — 1] ; 

f\ := BerlekampMassey(s); 

ifdeg(/i)=Dthen 

H := Hd{s) II Construct the Hankel matrix 

for i = 2, ...,n do 

b:=((r 7 ^-e):j = 0,...,D-l/; 
Compute c= (ci,...,c/))' :=H~ l b; 

Ji ■— Lk=0 c k+l x l> 

end 

return [f\,x 2 - f 2 ,- ■-,*„- fn}\ 
else 
[ return Fail; 
end 



3.1.2 Complexity 

In this complexity analysis and others to follow, we assume that the multiplication 
matrices are all known and neglect their construction cost. 

Suppose the number of nonzero entries in T\ is iVi. The Wiedemann algo- 
rithm (both construction of the linearly recurring sequence and computation of its 
minimal polynomial with the Berlekamp-Massey algorithm) will take 0{D (N\ + 
log(D))) field operations to return the minimal polynomial f\ |[36ll . 

Next we show how the linear system © can be generated for free. Note that for 
any a,b € K^ Dxl ^> and T E K^ DxD \ we have (a, Tb) = (T'a, b), where T' denotes 
the transpose of T. Thus in (f5]) and © 

(r,T{e) = ((r/)'r,e), <r,7>,) = {(T[) j r, Vi ). 

Therefore, when computing the sequence ©, we can record {T[)'r (i = 0, . . . , 2D — 
1) and use them for construction of the linear equation set ©. 

First, as each entry (r, T l +J e) of the Hankel matrix H can be extracted from 



10 



the sequence ©, the construction of H is free of operations. What is left now is 
the computation of ((Tf)Jr,Vi), where (T()->r has already been computed and Vj = 
Tje = NormalForm(x,). Without loss of generality, we can assume that NormalForm(x,) = 

x\ (this is not true only if there is a linear equation x,- H in the Grobner basis G\ , 

and in that case we can eliminate the variable xf). Consequently Vt is a vector with 
all its components equal to except for one component equal to 1. Hence comput- 
ing ((T[)ir,Vi) is equivalent to extracting some component from the vector {T[)ir 
and there is not additional cost. 

For each i = 2,...,n, solving the linear equation set He = 6, only needs 0(Dlog(D)) 
operations if fast polynomial multiplication is used [8]. Summarizing the analyses 
above, we have the following complexity result for this method. 

Theorem 3.2 Assume that T\ is constructed (note that Tz, ■ ■ ■ , T n are not needed). 
If the minimal polynomial of ® computed by the Berlekamp—Massey algorithm is 
of degree D, then the complexity of this method is bounded by 

0(D(Ni +log(D)) + (n- l)Dlog(D)) = 0(D{N x + nlog(D))). 



This complexity almost matches that of computing the minimal polynomial of 
the multiplication matrix T\ if n is small compared with D. 

3.1.3 Example 

We use the following small example to show how this method applies to ideals in 
shape position. Given the Grobner basis of a O-dimensional ideal / C Fn [x\ ,X2,x^\ 
w.r.t. DRL 

d = \x\ + 9x 2 + 2xi + 6, x\ + 2x 2 + 9, x 3 + 9] , 

we first compute the degree of / as D = 4, the canonical basis B = [l,xi ,X2,x\X2\, 
and the multiplication matrices T\, Tj and T^. 

With the random vector r = (8,4,8,6)' G K^ 4xl \ we can construct the linearly 
recurring sequence 

s= [8,4,0,7,6,8,10,10]. 

Then the Berlekamp-Massey algorithm is applied to s to obtain the minimal poly- 
nomial /i = x\ + 8xi + 9. From the equality deg(/i ) = D = 4, we know now the 
input ideal / is in shape position. 



The Hankel coefficient matrix 
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is directly derived from s. Next take the computation of the polynomial X2 — 
f2(xi) G G2 for example, the vector b = (8,6,8,3)' is constructed. The solution 
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of the linear equation set He = b being c = (1,0,5,0)', we obtain the polynomial 
in G 2 as X2 + 6x 2 { + 10. The other polynomial xj, — fy(xi) can be similarly com- 
puted. In the end, we have the Grobner basis of / w.r.t. LEX 

G 2 = [4 + 8xi + 9, x 2 + 6^ + 10, x 3 +9]. 

3.2 Deterministic algorithm to compute Grobner basis of radical of 
the ideal 

As already explained in Remarks |3. II the classical Wiedemann algorithm is a prob- 
abilistic one. For a vector chosen at random, it may only return a proper factor f\ 
of the polynomial f\, i.e., f\\f\ but f\ 7^/1. In fact, the deterministic Wiedemann 
algorithm can be applied to obtain the univariate polynomial f\ , then one knows 
for sure whether / is in shape position or not. The main difficulty is to compute the 
other polynomials / 2 , . . . ,f n in a deterministic way. 

In the following we present an algorithm to compute the Grobner basis of the 
radical of the ideal /. Indeed, in most applications, only the zeros of a polynomial 
system are of interest and we do not need to keep their multiplicities. Hence it is 
also important to design an efficient method to perform the change of ordering of 
Grobner basis of an ideal / in a way that the output is the Grobner basis of y/l. 

3.2.1 Deterministic version of the Wiedemann algorithm 

The way how this deterministic variant of the Wiedemann algorithm proceeds is 
first recalled. Instead of a randomly chosen vector in the classical Wiedemann 

algorithm, in the deterministic version all the vectors of the canonical basis of 

K (0xl) 

ei = (l,0,...,0y,e 2 = (0,l,0,...,0) f ,...,e D = (0,...,0,iy 

are used. One first computes the minimal polynomial f\ t \ of the linearly recurring 
sequence 

[( ei ,r/e):j = 0,...,2D-l]. (7) 

Suppose d\ = deg(/ii), and b\ = f\^(T\)e. If by = 0, one has f\\ = f\ and the 
algorithm ends; else it is not hard to see that the minimal polynomial f\^ of the 
sequence 

[(e 2 ,r/6i):; = 0,...,2(D-Ji)-l] 

is indeed a factor of /i//i,i, a polynomial of degree <D — d\ (that is why only the 
first 2{D — d\) terms are enough in the above sequence). Next, one computes 6 2 = 
/ii/i j2 (r)e and checks whether 6 2 = 0. If not, the above procedure is repeated 
and so on. This method ends with r (< D) rounds and one finds f\ = f\\ ■ ■ ■ fi >r . 
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3.2.2 Deterministic algorithm description 

First we study the general case when a factor of f\ is found. Suppose a vector 
w £ IC( Dxl ) is chosen to construct the linearly recurring sequence 

[(w,T{e}:i = 0,...,2D-l}, (8) 

and the minimal polynomial of this sequence is f\ , a proper factor of f\ of degree 
d. We show how to recover the Grobner basis of 7 + (f\ ) w.r.t. LEX. Since the ideal 
I is in shape position, it is not hard to see that the ideal I + (fi) is also in shape 
position, and its Grobner basis w.r.t. LEX is indeed [fi,X2 — fi, ■ ■ ■ ,x n — f n ], where 
fi is the remainder of /,• modulo f\ for i = 2, . . . , n. 

Now for each i, we can construct the linear system similar to © 

d-\ 
(w,TJT i e) = 1 £yk-{w,T 1 k+i e), j = 0,...,d-l, (9) 

where yg , . . . , yd- 1 are the unknowns. As the dxd Hankel matrix of © is invertible 
by Theorem 13. II there is a unique solution c r -,o,c,-i , . . . ,c r -^_i for ©. Next we will 
connect this solution and a polynomial in the Grobner basis of I + (f\), and the 
following lemma is useful to show this connection. 

Lemma 3.3 Suppose f\ is the minimal polynomial of ^for some w £ K(° x1 ), t\ 
the multiplication matrix ofx\ of the ideally (/i) w.r.t. <i, and e = (1,0,... ,0) £ 
j^(rfxi) tne canon i ca i basis of 1 in K[xi,...,x n ]/(I -\- (/i)). Then f\ is also the 
minimal polynomial of [e, T\e, T^e, . . .]. 

Proof Suppose fy = xf + Y?k=o akX \- Then according to FGLM criteria, for the 
ideal /+(/i), 

Tfe = £ a k t\e 

k=0 

is the first linear dependency of the vectors e,7ie,7f e,... when one checks the 
vector sequence [e, T\e, Tfe, ...]. That is to say, fi is also the minimal polynomial 
of [e,7\e, ife, . ..]. □ 



Proposition 3.4 Suppose w £ KS Dx l > is such a vector that a proper factor f\ of f 



of degree d < D is found from the linearly recurring sequence ([8]>. Then for each 
i = 2,...,n, the polynomial x,- — £fZn C[^x\, where c,-,o, c^\ , . . . , c^d-i is the unique 
solution of ©, is in the Grobner basis ofl+ {f\) w.r.t. LEX. 

Proof Let 7\ , . . . , 7^ be the multiplication matrices of the ideal I + (fi) w.r.t. < i . 

For each i = 2, . . . , n, suppose x,- — E*=n Cijtx\ is the corresponding polynomial 
in the Grobner basis of I + (J\ ) w.r.t. LEX. Then 7)e = YftZo ^yt^e holds, and for 
any vector w £ K^ x l \ we have 

d-\ 
(w,f/7-e) = Y, £ i,k-(w,f^ +J e), j = 0,...,d-\. 



k=0 
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As long as w is chosen such that the coefficient matrix is invertible, the coefficients 
Cifl > cv, 1 1 • • • > Ci,d- 1 will De the unique solution of the linear equation set 

d-\ 

(w, f(tie) = £ y k ■ (w,T, k+j e) , j = 0,...,d-\. (10) 

k=0 

Therefore, to prove the correctness of the proposition, it suffices to show that 
there exists w G K^ x l > such that the coefficient matrix of (fTOl is invertible, and 
that the two linear equation sets © and (flOl share the same solution. In particular, 
we will prove © and (flOl are the same themselves for some w. 

To prove that, we need to show the two Hankel matrices and the vectors in the 
left hands of © and (flOl are the same. That is, for some vector w 

(i) (w,T^e) = (w,T^e), for j = 0, ...,2d — 2; 

(ii) (w, r/zje) = (w, f(fie), for j = 0, . . . , d- 1. 

Next we will prove these two arguments respectively, 
(i) We take the first d equations in (i) 

{wj(e) = {w,f(e), j = 0,...,d-l. 

As the vectors e,7\e, .. . ,!Tj _1 e are linearly independent, the above linear equation 
set has a unique solution w for the unknown w. From Lemma [331 the vector 
sequence [e, 7\ e, 7\ e, . . .] and the sequence ([U) share the same minimal polynomial 
/i of degree d. Thus there exist ao> • • • j^d-i £ K such that 

rf-i rf-i 

f } d e = £ akTfe, (w, Tfe) = £ a k {wj?e). 

k=0 k=Q 

Hence 

d- 1 d- 1 rf- 1 

(w,ffe) = (w, £ a*?^) = £ a k {w,f^e) = £ a k {wrfe) = (w,Tfe). 

k=Q k=0 k=0 

Other equalities in (i) for j = d + 1 , ■ • • , 2d — 2 can also be proved similarly. Actu- 
ally, the equality (to, r/e) = {wq, f(e) holds for any j = 0, 1, 

(ii) Since there is a polynomial *,• — EfcTn a 'k x \ m tne Grobner basis of / w.r.t. 
LEX, where a f , . . . ,a' D _ l G 1C, we know 7}e = X^T a' k T^e. Then on one hand, for 
the vector w and any i = 0, . . . , d — 1 , we have 

<™,r/7}e>=£^<™,7f + ; e >. 

/t=0 

On the other hand, as x- { — Y^=o a 'k x \ ^ ^> we nave x < ~~ ££o a £-*i ^ ' + (/i ) > an( l 
thus 7je = Ylk=o a'kT\&- Therefore for the vector w and any j = 0, . . . , d — 1, 

(W,f/2jg) = ^a'^f^e) = Z<S k («>,T} +i e) = (rv,T^T t e). 
This ends the proof. □ 
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Now let us return to the special case of the deterministic Wiedemann algorithm, 
where unit vectors are used to find f\ = f\\ ■ •• f\. r with r <D and deg(/i,) = <f,. 
Suppose deg(/i) = D so the ideal / is verified in shape position. In the z'th step 
of the algorithm, the unit vector e, is applied to construct the linearly recurring 
sequence 



;"-l 



[{e i ,T{b i _ 1 ):j = 0,...,2(D-Y[d k )-l], 

k=\ 

where 6,-_i = Y\!iZ}\f\,k{T\)e. With this sequence the factor f\ j is computed. As 
the above sequence is the same as 

[(([\fi,k(Ti)ye i ,T^e):j = o,...,2(D-fld k )-i], 
k=\ k=\ 

from Proposition 13 .4 1 we can recover efficiently the Grobner basis of /+ (/i,j) w.r.t. 
LEX by constructing and solving linear equation sets with Hankel coefficient ma- 
trices. 

So we have at hands the factorization /i =/i 1 • • • fi, r , together with the Grobner 
basis of / + (fij) w.r.t. LEX for i = 1 , . . . , r. Suppose the Grobner basis for i is 

Pi = [f\,i,X2 ~ flh •■•,%- fn,i] ■ (1 1) 

Then to recover the polynomials fj in ((U) for j = 2, . . . ,n, we have the following 
modulo equation set constructed from P\,.. . ,P r : 



'fj-fj,l m °d/i,i 



Jj = fjs mod/i, r 



(12) 



Now it is natural to give a try of the Chinese Remainder Theorem (short as CRT 
hereafter). 

To use the CRT, we have to check first whether f\\ ,... ,/i >r are pairwise co- 
prime. One simple case is when f\ is squarefree, or in other words the input ideal / 
is radical itself. In that case, the direct application of CRT will lead to the Grobner 
basis G of I w.r.t. LEX, and the change of ordering ends. 

When the polynomial f\ is not squarefree, the CRT does not apply directly. 
In this case, the Grobner basis of yl w.r.t. LEX is our aim. Before the study on 
how to recover this Grobner basis, we first make clear how a polynomial set of 
form (O can be split to a series of polynomial sets with a certain zero relation 
according to some factorization of f\ . The following proposition is a direct result 
of [24, Proposition 5(i)], and it is actually a splitting technique commonly used in 
the theory of triangular sets [351. In what follows, Z(F) denotes the common zeros 
of a polynomial set F G K[xi , . . . ,x n ] in K , where K is the algebraic closure of K. 

Proposition 3.5 Let T C "K[x\ ,... ,x n ] be a polynomial set in the form 
[h(xi),X2-t2(xi),...,x n -t n (xi)}, 
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and t\ = fi i • • ■ h tT . For i= 1 , . . . , r, define 

1 [}) — \f\,ii^2 '2,*') • ■ • ;•"-« <«,<J j 

where tjj is the remainder of tj modulo t\i for j = 2, . . . ,n. Then we have the 
following zero relation 

r 

Z(T) = [jZ(T(i)). (13) 

Let /j be the squarefree part of f\. As each Pj in ffTTTt satisfies the form in 
Proposition 13.51 we can compute t new polynomial sets Pj whose univariate poly- 
nomials in x\ is fi j for j = 1, . . . , t, such that f x = Ytj=\ f \ /> an d /i / are pairwise 
coprime. These new polynomial sets can be found in the following way. Set p = f x . 
We start with j = 1 and computes /j ■ = gcd(/ij,p). As long as this polynomial 
is not equal to 1, a new polynomial set Pj whose univariate polynomial is f i ■ is 
constructed from Pj by Proposition 13.51 Next set p := p/f\ , and check whether 
p = 1. If so, we know we already have enough new polynomial sets; otherwise 
j := j + 1, and the process above is repeated. 

Now we reduce the current case to the earlier one with f x squarefree and 
P\,...,P t to construct the modulo equation sets. Thus the Grobner basis of \fl 
w.r.t. LEX can be obtained similarly (note that extracting the squarefree part of f\ 
results in the radical of /). 

The whole method based on the deterministic Wiedemann algorithm is summa- 
rized in Algorithm [3] below. The subfunction SqrfreeQ returns the squarefree part 
of the input polynomial. The operator "cat" means concatenating two sequences. 

Remark 3.2 If the factors /i,i, ■ • . ,/i,r of /i returned by the deterministic Wiede- 
mann algorithm are pairwise coprime (which needs extra computation to test), the 
Grobner basis of / w.r.t. LEX can be computed from the CRT. 

The method of the deterministic version described above is also applicable to 
the Wiedemann algorithm with several random vectors. To be precise, when the 
first random vector does not return the correct polynomial f\ , one may perform a 
similar procedure as the deterministic Wiedemann algorithm by updating the se- 
quence with a newly chosen random vector (instead of e,- in the basis) and repeat- 
ing [36]. In that case, the method above with CRT can also be used to compute the 
Grobner basis of a/7 w.r.t. LEX. 

3.2.3 Complexity 

Next the computational complexity, namely the number of field operations needed, 
for the deterministic method for ideals in shape position is analyzed. 
(1) In total the deterministic Wiedemann algorithm needs 

0(D(M+Dlog(D)loglog(D))) 
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Algorithm 3: Shape position (deterministic) G 2 '■= ShapeDet(Gi, <i) 

Input: Gi , Grobner basis of a O-dimensional ideal / C K.[xi, . . . ,x n ] w.r.t. < 1 
Output: G 2 , Grobner basis of \fl w.r.t. LEX if / is in shape position; Fail, 
otherwise. 

Compute the canonical basis of K[x\ ,... ,x n ]/(G\) and multiplication 

matrices T\,...,T n ; 

e 1 = (l,0,...,Oy, e2 = (0,l,0,...,0) f ,...,ez) = (0,...,0,iyGlK( Z5xl ); 

k:=l; F:=[}; /;= 1; d := 0; b = e i; S = []; 

while b / do 

s:=[{e k ,T{b) :i = 0,1,..., 2(n-d)-l}; 
g := BerlekampMassey(i'); 

f:=f-g; d:=deg(f); F:=Fcat [g]; 6:=g(7i)6; S := Scat [j]; 
Jt:=Jt+l; 
end 

(Suppose F = [/ M , . . . ,f Lr ]) /1 := nLi h,» 
ifdeg(/i)^Dthen 
I return Fail 
else 

for / = l,...,r do 
d,- :=deg(/ u ); 
for j = 2,... ,n do 

Construct the Hankel matrix // ; and the vector b from S; 
Compute c = (ci, . . . ,c rf ,) f := Hj l b\ fij := £fL Cfc + i4; 
end 
end 

J^Sqrfree^); 
if7i//ithen_ 

Compute { \f 1 j,x 2 - f 2 j, ■ ■ ■ ,*n - f n ,j] > J = 1, ■ • • >0 from 
{ l/i , j, *2 - /2,i , •_ L - ,*« - j»,i] = * = 1 > • • • > r } by Proposition [33] such 
that f 1 = Y['j=i fij an d fij are pairwise coprime; 
end 

for j =2,...,ndo 
I Solve the modulo equation set (fT2l) to get /_,•; 
end 

return [/j,^ - f%, . . . ,x„ - /„] 
end 



operations if fast polynomial multiplications are used [36]. Here N\ still denotes 
the number of nonzero entries in T\ . 

(2) Next at most D structured linear equation sets with Hankel coefficient ma- 
trices are constructed and solved, each with maximum operations 0(Dlog(D)). 
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Hence this procedure needs 0(D 2 log(D)) operations at most. 

(3) The squarefree part f l of f\ can be obtained with complexity at most 
0{D 2 log(D)) for the case when K has characteristic and 0{D 2 log(D) +D\og(q/p)) 
for characteristic p > respectively, where |K| = q 041 Theorem 14.20 and Exer- 
cise 14.30]. For the case when f\ is not squarefree, suppose t new polynomial sets 

P i , . . . , P t are needed, and deg (f u ) = d\ for i = 1 , . . . , t . To compute each set P, of 
the form d4), n — 1 polynomial divisions are needed to find the remainders, with 
complexity 0(ndiD). Hence the total complexity to obtain P\, . . . ,P t is 

0{j^ndiD) = 0{nDJ^di) < 0{nD 2 ), 
i=l j=l 

for we have Yli=\ di = deg (/J < D. 

(4) Solving the modulo equation set (fT2l for each j = 2, . . . ,n requires 0(D 2 ) 
operations at most by [34, Theorem 5.7]. Thus in total 0{nD 2 ) operations are 
needed for the CRT application. 

Therefore, we have the following complexity result for the method with the 
deterministic Wiedemann algorithm. 

Theorem 3.6 Assume that T\ is known. If the input ideal I is in shape position, 
then this deterministic method will return the Grobner basis of \JI w.r.t. LEX with 
the complexity 

0{D{N X +D(log(D)loglog(D) +«))). 

3.2.4 Example 

Here is a toy example to illustrate how the deterministic method works. Consider 
an ideal / in F2[jci ,x?\ generated by its Grobner basis w.r.t. DRL 

G\ := [x 2 x\+x\+xi + l,x\+x\+x 2 +'\-,x\+x 2 2 }. 

Its Grobner basis w.r.t. LEX is 

G 2 = [fi := (jci + lf{x\ +xi + l) 2 ,x 2 +x\ +x\ + 1], 

from which one can see that / is in shape position. 

From G\ the canonical basis B = [\,x\ ,x 2 ,x 2 ,xix 2 ,x i l ,x 2 x 2 ] and the multiplica- 
tion matrices T\ and T 2 are first computed. With a vector r = (1, 1,0, 1,0, l,0) f 6 
¥ 2 x generated at random, the classical Wiedemann algorithm will only return 
a proper factor (xi + l)(x 2 +x\ + 1) of f\, and whether / is in shape position is 
unknown. 

Next we use the deterministic Wiedemann algorithm to recover f\. With e\ = 
(1,0, . . . ,0)', a factor /u = (xi + \) 2 {x\ +x\ + \) of/i is found with the Berlekamp- 
Massey applied to the sequence ©. Then we update the vector 

6 = / u (7i)e = (0,1,1,0,0,0,0)', 
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and execute the second round with e2 = (0, 1,0, ... ,0)', obtaining another factor 
/i,2 = (xi + !)(*! +*i + 1). This time the updated vector 6 = 0, thus the determin- 



istic Wiedemann algorithm ends, and f\ is computed as /i, 1/1,2- As deg(/i ,1/1,2) = 
D, now / is verified to be in shape position. 

Then we construct the linear equation sets similar to © to recover /^i and /2,2 
respectively. The first one, for example, is 
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After solving them, we have the Grobner bases of /+ (/1 j) and /+ (/i,2) respec- 
tively as 

Pi = [(xi + l) 2 (xj +Xi + l),x 2 +xi], 
Pi = [(xi + \){x\ + X\ + 1),JC 2 +Xl]. 

Then the squarefree part f l of f\ is computed, and we find that / is not radical, 
and thus only the Grobner basis G2 of a/7 w.r.t. LEX may be computed. From 
/i,2 = /1. we directly have G% =Pz, and the algorithm ends. 

The way to compute G2 by CRT, which is more general, is also shown in the 
following. Two new polynomial sets 

Pi = [Xl + 1,X 2 + 1], P% = [*i+*i + 1,^2+Xi] 



are first computed and selected according to f 1 by Proposition 13.51 Then the mod- 
ulo equation set 

{fi = x 2 + 1 m °d x \ + 1 , 
fi = x 2 + x \ m °d x\+x\ + \ 

as (TT2T > is solved with CRT, resulting in the same G2. One can check that G2 is the 
Grobner basis of yl w.r.t. LEX with any computer algebra system. 

3.3 Incremental algorithm to compute the univariate polynomial 

For a 0-dimensional ideal / C K[jci, . . . ,x n ], the univariate polynomial in its Grobner 
basis w.r.t. LEX is of special importance. For instance, it may be the only poly- 
nomial needed to solve some practical problems. Furthermore, in the case when 
K is a finite field, after the univariate polynomial is obtained, it will not be hard 
to compute all its roots, for one can simplify the original polynomial system by 
substituting the roots back, and sometimes the new system will become quite easy 
to solve. 

Besides the two methods in the previous parts, next the well-known incremen- 
tal Wiedemann algorithm dedicated to computation of the univariate polynomial is 
briefly recalled and discussed. 
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In the Wiedemann algorithm, the dominant part of its complexity comes from 
construction of the linearly recurring sequence (0(DN\)), while the complexity 
of the Berlekamp-Massey algorithm is relatively low (0(Dlog(D))). Hence the 
idea of the incremental method is to construct the sequence incrementally to save 
computation and apply the Berlekamp-Massey algorithm to each incremental step. 

We start with the linearly recurring sequence [(r,T{e) : i = 0, 1] and compute 
its minimal polynomial with the Berlekamp-Massey algorithm. Next we proceed 
step by step with the sequence 

[(r,r/e):» = 0,...,2*-l] 

until the returned polynomial coincides with the one in the previous step. Then this 
minimal polynomial equals the univariate polynomial / we want to compute with 
a large probability. 

Suppose deg(/) = d. Then the number of steps the method takes is bounded 
by d + 1. In other words, the method stops at most after the sequence [(r, T{e) : i = 
0, . . . ,2d+l] is handled. The number of field operations to construct the sequences 
is 0(dN[), while the total complexity to compute the minimal polynomials with the 
Berlekamp-Massey algorithm is Ofck={ ^ 2 ) = 0(d 3 ) (note that in the incremental 
case, the fast Berlekamp-Massey with complexity 0(klog(k)) is not applicable). 
Therefore the overall complexity for the incremental Wiedemann method to com- 
pute the univariate polynomial is 0{dN\ +d 3 ). As can be seen here from this 
complexity, this incremental method is sensitive to the output polynomial /. When 
the degree d is relatively small compared with D, this method will be useful. 

4 General case: B MS -based algorithm 

In the general case when the ideal / may not be in shape position, perhaps those 
methods described in Section [3] will not be applicable. However, we still want to 
follow the idea of constructing linearly recurring sequences and computing their 
minimal polynomials with the Berlekamp-Massey algorithm. The way to do so is 
to generalize the linearly recurring sequence to a multi-dimensional linearly recur- 
ring relation and apply the BMS algorithm to find its minimal generating set. 

4.1 Algorithm description 

We first define a ^-dimensional mapping E : Z> — > K as 

(s u ..., Sn )^{r,T?»-T*>e), (14) 

where r € KA Dxl ) is a random vector. One can easily see that such a mapping is a 
^-dimensional generalization of the linearly recurring sequence constructed in the 
Wiedemann algorithm. 

Note that T* 1 ■ ■ ■ T*"e in the definition of E above is the coordinate vector of 
(s\,...,s„) in the FGLM algorithm. As a polynomial / in the Grobner basis of / 
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is of form ©, and the linear dependency CQ holds, one can verify that / satisfies 
© and thus is a polynomial in 1(E). The BMS algorithm is precisely the one 
to compute the Grobner basis of 1(E) w.r.t. to a term ordering, so one may first 
construct the mapping E via T\ , . . . , T„, and attempts to compute the Grobner basis 
of / from the BMS algorithm applied to 1(E). 

We remark that / is in 1(E) for any vector r. In fact, the idea above is a multi- 
dimensional generalization of the Wiedemann algorithm. The minimal polynomial 
g of the Krylov sequence [b,Ab,A 2 b, . . .] is what the Wiedemann algorithm seeks, 
for g directly leads to a solution of the linear equation Ax = b for a nonsingular 
matrix A and vector b. Then a random vector is chosen to convert the sequence to 
a scalar one 

[(r,b),(r,Ab),(r,A 2 b),...}, 

and the Berlekmap-Massey algorithm is applied to find the minimal polynomial of 
this new sequence, in the hope that g can be obtained. While the method proposed 
here converts the mapping from (s\, . . . , s n ) to its coordinate vector in the FGLM to 
a ^-dimensional scalar mapping with a random vector, and then the BMS algorithm 
(generalization of Berlekamp-Massey) is applied to find the minimal polynomial 
set, which is also the Grobner basis, w.r.t. to a term ordering. 

This method for computing the Grobner basis of / makes full use of the spar- 
si ty of T\ , . . . ,T n , in the same way as how the Wiedemann algorithm takes advan- 
tage of the sparsity of A. The method is a probabilistic one, also the same as 
the Wiedemann algorithm. This is reasonable for the ideal 1(E) derived from the 
^-dimensional mapping may lose information of / because of the random vector, 
with / C 1(E). Clearly, when / is maximal (corresponding to the case when g in 
the Wiedemann algorithm is irreducible), 1(E) will be equal to /. Furthermore, 
as polynomials in the Grobner basis are characterized by the linear dependency in 
CQ), we are always able to check whether the Grobner basis of 1(E) returned by the 
BMS algorithm is that of /. 

Remark 4.1 When the term ordering in the BMS algorithm is LEX, computation 
of the univariate polynomial in this method is exactly the same as that described in 
Section |3~T1 This is true because for the LEX ordering (xi < --• < x n ), the terms 
are ordered as 

[l,JCl,xf,... ,X2,X[X2,X l X2,- ■■], 

hence the first part of E is E((p\,0, . .. ,0)) = (r,Tf'e), and the BMS algorithm 
degenerates to the Berlekamp-Massey one. 

Another fact we would like to mention is that the BMS algorithm from Coding 
Theory is mainly designed for graded term orderings like DRL, for such orderings 
are Archimedean and have good properties to use in algebraic decoding [13 1. But 
it also works for other orderings, though extra techniques not contained in the 
original literature have to be introduced for orderings dependent on LEX (like LEX 
itself or block orderings which break ties with LEX). 
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Take the term ordering LEX for instance, an extra polynomial reduction is 
performed after every BMSUpdate() step to control the size of intermediate poly- 
nomials. This is actually not a problem for orderings like DRL, for in that case 
the leading term of a polynomial will give a bound on the size of terms in that 
polynomial. We also have to add an extra termination check for each variable xu 
otherwise the BMS algorithm will endlessly follow a certain part of the terms. For 
example, all variables in the sequence [l,xi ,JCj , . . .] are smaller than X2, and the 
original BMS does not stop handling that infinite sequence by itself. 

With all the discussions, the algorithm is formulated as follows. The "Termi- 
nation Criteria" here in this description mean that F does not change for a certain 
number of iterations. The subf unction Reduce (F) performs reduction on F so that 
every polynomial / G F is reduced w.r.t. F \ {/}, and IsGB(F) returns true if F is 
the Grobner basis of / w.r.t. LEX and false otherwise. 

Algorithm 4: General case G 2 := BMSbased(Gi, <i) 
Input: Gi , Grobner basis of a O-dimensional ideal / C K[jci , . . . ,x n ] w.r.t. < i 
Output: Grobner basis of/ w.r.t. < 2 ; or Fail, if the BMS algorithm fails 
returning the correct Grobner basis 

Compute the canonical basis of K[xi, . . . ,x n ]/(G\) and multiplication 

matrices T\,...,T n ; 

Choose r G K^ 5 * 1 ) at random; 

u:=0; F:=[l]; G:=[]; E:=[]; 

repeat 

e:=(r,T l Ul ---T^e); 

E := E cat [e\\ 

F,G:= BMSUpdate(F,G,u,£); 

u := Next(tt) w.r.t. < 2 ; 

F := Reduce (,F); 
until Termination Criteria; 
ifIsGB(F)then 
I return F 
else 
I return Fail 
end 

The correctness of Algorithm |4] is obvious. Next we prove its termination. 
Once the loop ends, the algorithm almost finishes. Hence we shall prove the termi- 
nation of this loop. Clearly when the polynomial set F the BMS algorithm main- 
tains turns to the Grobner basis of 1(E) w.r.t. <% the current termination criterion, 
namely F keeps unchanged for a certain number of passes, will be satisfied. And a 
sufficient condition for F being the Grobner basis is given as Theorem !4.4l below. 
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4.2 Complexity 

Part of earlier computation of values of E can be recorded to simplify the com- 
putation at Line|4] Suppose the value of E at a certain term (u\, «2, • • ■ ; «/-i ; Ui — 
l,Ui+\, . . . ,u n ) 

e = T^---T™'- l ---T™»e 

has been computed and recorded. Then we know the value at u = (wi, ■■■,u n ) is 

<r,C...C»e> = (r,7}e>, 

for all 7] and Tj commute. Thus the computation of one value of E can be achieved 
within 0(N) operations, where N is the maximal number of nonzero entries in 
matrices T\,...,T n , 

Next we focus on the case when the target term ordering is LEX. The complex- 
ities of the three steps in Algorithm [T] are analyzed below. 

(1) As an extra reduction step is applied after each iteration, the numbers of 
terms of polynomials in F are bounded by D + 1. Denote by N the number of 
polynomials in G2. Then checking whether F is valid up to Next(iz) needs 0{ND) 
operations. 

(2) The computation of the new delta set A(Next(u)) only involves integer 
computations, and thus no field operation is needed. 

(3) Constructing the new polynomial set F + valid up to Next(tt) requires 
0{ND) operations at most. The readers may refer to iBTl H3l for the way to con- 
struct new polynomials. 

In step (1) above, new values of E other than e may be needed for the ver- 
ification. The complexity for computing them is still 0(N), and this is another 
difference from the original BMS algorithm for graded term orderings. After the 
update is complete, a polynomial reduction is applied to F to control the size of ev- 
ery polynomial. This requires 0(NND) operations, where ,/V denotes the maximum 
term number of polynomials in G^- To summarize, the total operations needed in 
each pass of the main loop in Algorithm |4]is 

0(N + ND + NND) = 0{N + NND) . 

Hence to estimate the whole complexity of the method, we only need an upper 
bound for the number of passes it takes in the main loop. 

Theorem 4.1 Suppose that the input ideal I C M.[x\ , . . . ,x n ] is of degree D. Then 
the number of passes of the loop in Algorithm^is bounded by InD. 

Before giving the proof, we need to introduce some of the proven results on the 
BMS algorithm for preparations. Refer to flT] [13] for more details. 

Denote the previous term of u w.r.t. < by Pre(-u). Given a ^-dimensional array 
E, suppose now a polynomial / € K[;q ,... ,x n ] is valid for E up to Pre(it) but not 
to u. Then the term u — lt(/) is called the span of / and denoted by Span(/), 
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while the term u is called the. fail of / and written as Fail(/). When / G 1(E), f 
is valid up to every term, and in this case we define Span(/) := oo. The following 
proposition reveals the importance of spans. 

Proposition 4.2 ( Corollary 9]) A(E) = {Span(/) : / 1(E)}. 

Define I(u) := {/ G K[;ci,. . . ,x n ] : Fail(/) > u}. Such a set is not an ideal 
but is closed under monomial multiplication: supposing that F G (a)(u), we have 
tF G (a) (u) for every term t G K[xi , . . . ,x n ]. 

Proposition 4.3 ([7, Proposition 6]) For each u, A(w) = {Span(/) : / I(u)}. 
Furthermore, v G A(tt) \ A(Pre(it)) if and only if v -< u and u — v G A(it) \ 
A(Pre(«)). 

The above proposition states when a term in A (is ) is determined, and it is going 
to be used extensively in the sequel. Also from this proposition, one can derive the 
following termination criteria for the BMS algorithm, which are mainly designed 
for graded term orderings like DRL. 

Theorem 4.4 ( [13, pp.529, Proposition (3.12)] ) Let c max be the largest element 
ofA(E) and s max be the largest element of{lt(g) : g G G}, where G is the Grobner 
basis of 1(E) w.r.t. <. 

(1) For all u > c max + c max , A(w) = A (is) holds. 

(2) For all v > c max + max.{c max ,s max }, the polynomial set F the BMS algorithm 
maintains equals G. 

As explained in Section 14.11 actually the term ordering LEX is not the one 
of interest in Coding Theory and does not possess some properties needed for a 
good order domain. But the results stated above are still correct. In particular, 
Theorem 14.41 indicates when the iteration in the BMS algorithm ends. For graded 
term orderings like DRL, once the termination term is fixed, the whole intermediate 
procedure in the BMS algorithm is also determined. However, for LEX it is not the 
case. We have to study carefully what happens between the starting term 1 and the 
termination term indicated by Theorem 14.41 

Next we first illustrate the procedure for a 2-dimensional example derived from 
Cyclic5. Both the delta set (marked with crosses) and the terms handled by the 
BMS algorithm (with diamonds) are shown in Figure Q] 

The Cmax and s max in Theorem 14.41 are respectively (4,6) and (0,7). In fact, 
the BMS algorithm obtains the whole delta set at (8, 12) = c max + c max , and the 
polynomial set it maintains grows to the Grobner basis at (4, 13) = c max + s max , 
which is also where the algorithm ends. 

Next we go into some details of what happens when a diamond row is handled 
by the BMS algorithm. We call a diamond (or cross) row the jth diamond (or cross) 
row if terms in this row are (i, j). Then for the 0th diamond row, the BMS algorithm 
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Figure 1: Delta set (+) and terms needed (o) for Cyclic5-2 

degenerates to the Berlekamp-Massey one to compute the univariate polynomial 
f\(x\). Here 30 diamond terms are needed because the minimal polynomial is of 
degree 15. 

For other rows in Figure [Q from Proposition 14-3 J one knows that at a 7th di- 
amond rows with an odd j, the delta set does not change. Thus such diamond 
rows are only bounded by the latest verified row in the delta set. This is because 
otherwise a wrong term in the delta set will be added if other diamond terms are 
handled. For example, the 3rd diamond row is of the same length as the 1st cross 
row, while the 5th diamond row is as that of the 2nd cross one. 

For a 2kth diamond row, its number is related to two cirteria. On one hand, 
again from Proposition 14. 3 1 the kth cross row is determined while the 2/cfh diamond 
row is handled in the BMS algorithm. Denote by c max {k) the largest term in the 
kth cross row, then terms up to c max {k) +c max (k) in the 2kth diamond row have to 
be handled to furnish the kth cross row. On the other hand, the number of 2kth 
diamond row is also bounded by the latest verified cross row, as the odd diamond 
ones. The first criterion is shown by the 6th diamond and the 3rd cross rows, while 
the 4th diamond row is the result of both criteria. 

For a term u = (u\, . .. ,k,-,0, . . . ,0) 6 K[xi,.. . ,xi\, in the proof below we write 
it as u = (u\,. . . ,ut) for simplicity, ignoring the last n — i zero components in the 
terms. 

Proof (of Theorem 14. II ) Suppose G is the Grobner basis of 1(E) the BMS algo- 
rithm computes. Denote the number of terms needed in the BMS algorithm to 
compute G n K[x\ ,... ,jc,-] by %i, and A,- := A (is) n K[xy ,... ,Xj]. From I C 1(E) 
one knows that A (is) is a subset of the canonical basis of K[x\,...,x n ]/I, thus 
|A(is)| < D. Therefore to prove the theorem, it suffices to show 2n|A(is)| is an 
upper bound. 

We induce on the number of variable / of K[x\ ,... ,x,-]. For i = 1, the BMS 
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algorithm degenerates to the Berlekamp-Massey, and one can easily see that %\ < 
2|Ai| holds. Now suppose %k < 2&|Afc| for k(< n). Next we prove %k+\ < 2(k-\- 
l)|A fe+ i|. 

As previously explained, in the terms to compute GflKfxi, . . . ,X£ + i], the terms 
(wi,... ,Uk,2l) are determined by two factors: terms (vi,...,Vfc,/) in Ajt+i, and the 
latest verified terms in the delta set. First we ignore those (ui,...,Uk,2l) terms 
determined by the latter criterion, and denote by .5,1+1 all the remaining ones in 
K[xi,. .. ,J<kfi]. We claim that |^+i| is bounded by (2k + l)|Afc+i|. 

From Theorem 14.41 we can suppose there exists some integer m, such that 

3Tk+\= (J $k+i,u A fe+ i= |J Afc+ij, (15) 

Z=0,...,2m+1 7=0,.. .,m 

where 

5^+1,/ :={«£ Jfc+i :u = (ui,...,Uk,l)}, 
Ak+ij :={«G Ajt+i :« = (Mi,...,M t ,7)}. 

Clearly |^+i,o| = Xk < 2k\Aj c \, and A^+i.o = A/;. One can see that \^k+l,2j\ 
is bounded by either 2&|Afc| = 2/clA^+^ol (if j = 0) or 2|A^ + ij| (< 2k\Ai c+ ij\). 
Furthermore, \^l+i,2j+i\ * s bounded by |Ajt+ij|, the number of the latest verified 
delta set. Hence we have 

\&k+i,2j\ + |^+i,2y+i | < (2k+l)\Ak+ij\, 

which leads to |«5£+i| < (2JtH-l)|A fc+ i|. 

Now we only need to show the number of all the previously ignored terms, 
denoted by =^' +1 , is bounded by |At+i|. Suppose £?L\ = Uies ^k+l i» wnere 5 is a 
set of indexes with \S\ < m in (fT5T >. and 

^+y:={«e ^'+i :u=(ui,...,u k ,i)}. 

Then for each z, \&£< j ,-| is bounded by the number of the latest verified delta set, 
say |A^ + i p . |. Thus the conclusion can be proved if one notices [Jies^-k+i.pi Q 

A k+l . " ' n 

Theorem 4.5 Assume that Ti,...,T n are constructed. The complexity for Algo- 
rithm^ to complete the change of ordering is bounded by 0(nD(N + NND)), 
where N is the maximal number of nonzero entries in the multiplication matri- 
ces T\,...,T n , and N and N are respectively the number of polynomials and the 
maximal term number of all polynomials in the resulting Grobner basis. 

4.3 Example 

Consider the ideal / C F65521 [^1,^2] defined by its DRL Grobner basis (x\ < X2) 

Gi ={x2 + 2x]x 2 + 2lxl + IIX1X2 +4x^2 + 22x^ +9x2 + 17xix 2 + 19x?+ 

2x 2 + 19xi + 5,XiX2 + lOx^ + 12xjx 2 + 20xi +21, x\ + 15xf + 19xi + 3}. 
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Here F65521 [x\,X2\/ {G\) is of dimension 12. Its basis, and further the multiplication 
matrices T\ and T2, can be computed accordingly. 

Now we want to compute the Grobner basis G2 of / w.r.t. LEX. With a vector 

r = (6757,43420,39830,45356,52762, 17712, 

27676, 17194, 138,48036, 12649, 11037)' G F^? 

generated at random, the 2-dimensional mapping E is constructed. Then BMSUpdate() 
is applied term by term according to the LEX ordering, with A(u) and the poly- 
nomial set F valid up to u shown in Tabled] For example, at the term (4,0), the 
polynomial xf + 62681;q +41493 € F is not valid up to (4,0). Then the delta 
set is updated as {(0,0), (1,0), (2,0)}, and F is reconstructed such that the new 
polynomial x\ + 6268 l;cf + 35812*1 + 18557 is valid up to (4,0). 
The first polynomial in G2'. 

gi=xf + 15x? + 19;ti+3 

is obtained at the term (7,0). Next BMSUpdate() is executed to compute other 
members of 1(E) according to the remaining term sequence [x2 , x\ xz , . . . , x\ , x\x\ ,...,], 
until the other polynomial in G2'. 

g2 =x\ + lx]x 2 2 + \5x\x2 + 2x\ +9 

is obtained at (3,5). Now the main loop of Algorithm |4] ends. Then one can easily 
verify that {g u g 2 } C G 2 anddim(F 6 552i[*2,*i]/(gi,g2)) = 12, thus G 2 = {gi,gz}- 



Term u 


A(«) 


F: polynomial set valid up to u 


(0,0) 


(0,0) 


X\,X2 


(1,0) 





xi+ 65437, x 2 


(2,0) 


(0,0), (1,0) 


x\ + 65437xi + 21672, .v 2 


(3,0) 


— 


xf + 62681jc, +41493,x 2 


(4,0) 


(0,0), (1,0), (2,0) 


,^ + 62681xf + 35812x, + 18557, x 2 


(5,0) 





x\ +30688^ + 45566x1 +54643,x 2 


(6,0) 


(0,0), (1,0), (2,0), (3,0) 


x\ + 30688*3 + 20026xf + 45766x! + 5434,x 2 


(7,0) 





«1>X2 


(0,1) 





gi ,x 2 + 65034x3 + 24330xf + 14876xi + 52361 


(1,1) 





gi ,x 2 + 64550x5 + 37707 xf + 48745*i + 7628 


(2,1) 





gl ,x 2 + 38842x5 + 5603x2+45755x! +44311 


(3,1) 





gi ,x 2 + 9449x3 + 20826x2 + 39078x! + 38885 


(0,2) 


(0,0), (1,0), (2,0), (3,0), (0,1) 


gi,x\ + 38885x 2 + 65360x5 + ln2x \ + 36000xi + 39469 
x 2 xi + 20826x] + 28385xf + 55917x, + 37174 



Table 1: Example for the BMS-based method 

Here is an example where this method fails. Let G = {x\,x\x2, x\x\,xQ C 
F65521 [jci,JC2]< Then the ideal (G) is 0-dimensional with degree D = 6. It is easy to 
see that G is Grobner basis w.r.t. both DRL and LEX. Starting from G as a Grobner 
basis w.r.t. DRL, the method based on the BMS algorithm to compute the Grobner 
basis w.r.t. LEX will not be able to return the correct Grobner basis, even the base 
field itself is quite large and different random vectors r are tried. 
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5 Putting all methods together: top-level algorithm 

In this section, we combine the algorithms presented in the previous parts of this 
paper as the following integrated top-level algorithm, which performs the change 
of ordering of Grobner bases to LEX. 

Algorithm 5: Top-level algorithm G 2 '■= TopLevel(Gi, <i) 
Input: G\ , Grobner basis of a O-dimensional ideal / C K[x\ ,...,x n ] w.r.t. < 1 
Output: G2, Grobner basis of / or yl w.r.t. LEX 
G 2 :=ShapePro(Gi,<i); 
ifG 2 / Fail then 
I return G 2 
else 

G 2 :=ShapeDet(Gi,<i); 
ifG 2 / Fail then 
I return G 2 
else 

G 2 :=BMSbased(Gi,<i); 
if G 2 / Fail then 
I return G 2 

else 
I return FGLM(G 1 ,< 1 ) 
end 
end 
end 

We would like to mention that to integrate these three algorithms, one needs 
to skip some overlapped steps in the three algorithms, like computation of the 
canonical basis and the multiplication matrices, and the choice of random vectors, 
etc. If one does not seek for the Grobner basis of \/7, that is to say, the multiplicities 
of the zeros are needed, then the deterministic invariant should be omitted. 

Thanks to the feasibility in each algorithm to test whether the computed poly- 
nomial set is the correct Grobner basis, this top-level algorithm will automatically 
select which algorithm to use according to the input, until the original FGLM one is 
called if all these algorithms fail. It is also a deterministic algorithm, though both 
the Wiedemann algorithm and the BMS-based method will introduce randomness 
and probabilistic behaviors to the individual algorithms. 

For an ideal in shape position, the probability for Algorithm [2] to compute the 
correct Grobner basis is the same as that of computing the correct minimal poly- 
nomial in the Wiedemann algorithm for one choice of a random vector, which has 
been analyzed in ll36l . When Algorithm [2] fails, the one based on the determinis- 
tic Wiedemann algorithm can tell us for sure whether the input ideal is in shape 
position, and return the Grobner basis of \fl. However, the probability for the 
BMS-based method to return the correct Grobner basis is still unknown. 
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6 Multiplication matrix T\ : sparsity and complexity 

In the previous description and complexity analyses of all the algorithms, the mul- 
tiplication matrices T\ , . . . , T„ are assumed known. In this section, for generic 
polynomial systems and the term ordering DRL, the multiplication matrix T\ is 
exploited, on its sparsity and cost for construction. We are able to give an ex- 
plicit formula to compute the number of dense columns in 7\, and we also analyze 
the asymptotic behavior of this number, which further leads to a finer complex- 
ity analysis for the change of ordering for generic systems. The term ordering is 
preassigned as DRL in this section without further notification. 

6.1 Construction of multiplication matrices 

Given the Grobner basis G of a O-dimensional ideal / w.r.t. DRL, let B = [ei , . . . , Er>] 
be the canonical basis of K[xi,. . . ,x„]/(G), and L := {lt(g) : g € G}. The three 
cases of the multiplication EiXj for the construction of the fth column of Tj in FGLM 
are reviewed below Ifl9l. 



(1) The term EjXj is in B: the coordinate vector of NormalForm(£,xy) is (0, . . . , 0, 1 , 0, . . . , 0) f , 
where the position of 1 is the same as that of epCj in B; 

(2) The term epcy is in L: the coordinate vector can be obtained easily from the 
polynomial g € G such that lt(g) = Sixy, 

(3) Otherwise: the normal form of EjXj w.r.t. G has to be computed to get the 
coordinate vector. 

Obviously, the /th column of Tj is sparse if case (a) occurs, thus a dense column 
can only come from cases (b) and (c). Furthermore, the construction for a column 
will not be free of arithmetic operations only if that column belongs to case (c). 
As a result, we are able to connect the cost for construction of the multiplication 
matrices with the numbers of dense columns in them. 

Proposition 6.1 Denote by M, the number of dense columns in the multiplication 
matrix Tj. Then the matrices T\,...,T„ can be computed within 0(D 2 ^" =1 M r ) 
arithmetic operations. 

Proof Direct result from the proof of Proposition 3. 1 in ITT91 . □ 

As shown in Section [3l among all multiplication matrices, T\ is the most im- 
portant one, and it is also of our main interest. However, for an arbitrary ideal, now 
we are not able to analyze the cost to construct T\ by isolating it from the others in 
Proposition 16.11 for the analysis on T\ needs information from the other matrices 
too. 

In the following parts we first focus on generic sequences which impose stronger 
conditions on 7\ so that the analyses on it become feasible. We show that the con- 
struction of T\ for generic sequences is free and present finer complexity results 
based on an asymptotic analysis. 
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6.2 Generic sequences and Moreno-Socias conjecture 

LetP= [/i,... ,f n ] be a sequence of polynomials inK[jfi, • • • ,x n ] of degree d\ , . .. ,d n . 
If d\ = ■ ■ ■ = d n = d, we call it a sequence of degree d. We are interested in the 
properties of the multiplication matrices for the ideal generated by P if /i, . . . ,/„ 
are chosen "at random". Such properties can be regarded generic in all sequences. 
More precisely, let U be the set of all sequences of n polynomials of degree d\,...,d n , 
viewed as an affine space with the coefficients of the polynomials in the sequences 
as the coordinates. Then a property of such sequences is generic if it holds on a 
Zariski-open in U. Next for simplicity, we will say some property holds "for a 
generic sequence" if it is a generic one, and also P is "a generic sequence" if its 
properties of our interest are generic. 

For a generic sequence [/i, ...,/„], its properties concerning the Grobner basis 
computation, in particular the canonical basis, are the same as [fi , . . . ,/*], where 
jf is the homogeneous part of /,• of the highest degree. That is to say, we only 
need to study homogeneous generic sequences, which are also those studied in the 
literature. Hence in the following part of this section, a generic sequence is further 
assumed homogeneous. 

Since we restrict to the situation where the number of polynomials is equal 
to that of variables, a generic sequence is a regular one (23). We first recall the 
well-known characterization of a regular sequence via its Hilbert series. 

Theorem 6.2 Let \f\ , . . . , f r ] be a sequence in K[x^ ,... ,x n ] with deg(f) = d{. Then 
it is regular if and only if its Hilbert series is 

nLi(i-^') 

(l-z)" 
Let P be a generic sequence of degree d. Then we know its Hilbert series is 

H(n,d):=(\-z d ) n /(l-z) n = (l+z + z 2 + ---+z d - [ )'\ (16) 

from which one can easily derive that the degree of (P) is d", and that the greatest 
total degree of terms in the canonical basis is (d — \)n. 

Grobner bases of generic sequences w.r.t. DRL have been studied in [27). A 
term ideal J is said to be a weakly reverse lexicographic ideal if the following 
condition holds: if t G J is a minimal generator of J, then J contains every term of 
the same total degree as t which is greater than t w.r.t. some term ordering. For the 
term ordering DRL, we have the following conjecture due to Moreno-Socias. 

Moreno-Socias conjecture ( |[26l0 Let K be an infinite field, and P = [f\ , . . . ,f n ] a 
generic sequence in K[xi,. . . ,x„] with deg(/,) = dj. Then \i{{P)), the leading term 
ideal of (P) w.r.t. DRL, is a weakly reverse lexicographical ideal. 

The Moreno-Socias conjecture is proven true for the codimension 2 case and 
for some special ideals for the codimension 3 case |[T][TT]. It has been proven that 
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this conjecture implies the Frtiberg conjecture on the Hilbert series of a generic 
sequence, which is well-known and widely acknowledged true in the symbolic 
computation community EHl . 

Proposition 6.3 Use the same notations as those in the Moreno-Socias conjecture. 
If this conjecture holds, then for a term u € lt((P)), any term v such that deg(ii) = 
deg(i>) and v > u is also in lt((P)). 

Proof If u is a minimal generator of lt((P)), then the conclusion is a direct result 
from the Moreno-Socias conjecture. Else there exists one minimal generator u^u 
such that u >~ u. As for any w such that deg(w) = deg(w) and w > u, we know 
w £ \t((P}). Then we can always find a term v G lt((P)) such that v y v. For 
example, construct v = u — u + v. If v is a term, then it suffices; otherwise the 
biggest term v such that deg(U) = deg(w) will work. This ends the proof. □ 

As Proposition 16.31 implies, the Moreno-Socias conjecture imposes a stronger 
requirement on the structure of the terms in lt((P)) for a generic sequence P. For 
the bivariate case, once a term u is known to be an element in lt((P)), the terms 
in lt((P)) determined by it are illustrated in Figure |2](left), and furthermore, in the 
right figure the shape all terms in lt((P)) form. 
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Figure 2: One term u E lt((P)) (•) and the terms it determines (■) / Minimal 
generators of lt((P)) (•) and terms in lt((P)) (■) 

The base field in the Moreno-Socias conjecture is restricted infinite. According 
to our preliminary experiments on randomly generated sequences over fields of 
large cardinality, we find no counterexample of this conjecture. As a result, we 
will consider it true and use it directly. The following variant of Moreno-Socias 
conjecture, which is more convenient in our setting, can be derived easily from 
Proposition [ 



Variant of Moreno-Socias conjecture Let K be an infinite field, P € K[x\ , 
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a generic sequence of degree d, and B the canonical basis o/K[xi,...,x„]/(P) 
w.r.t. DRL. Denote by B{k) the set of terms of total degree k in B. Then for 
k = 1, . . . , (d — \)n, B{k) consists of the first \B(k)\ smallest terms in all terms 
of total degree k. 



6.3 Sparsity and construction 

Let P C K[x\ ,... ,x n ] be a generic sequence, and G the Grobner basis of (P). Then 
polynomials in G can be assumed dense (in the case when IK is of characteristic 
or of large cardinality). As the number of dense columns in T\ will directly lead to 
a bound on the number of nonzero entries in T\ , the study of T\ sparsity is reduced 
to that of how many cases of (2) and (3) happen. Combining the Hilbert series of a 
generic sequence and our variant of Moreno-Socias conjecture, we are able to give 
the counting of the dense columns in T\ . 

Proposition 6.4 Let K, P, B and B(k) be the same as those in the Moreno-Socias 
conjecture variant. If the Moreno-Socias conjecture holds, then the number of 
dense columns in the multiplication matrix T\ is equal to the greatest coefficient in 
the expansion of(l+z-\ Vz^^Y- 

Proof Let k' = (d — \)n, and denote by 3f( k > be set of all terms in K[*i , . . . ,x n ] of 
total degree k. 

Suppose that u is the /th smallest term in ST^, Then^itt is still the /th smallest 
term in £?( k+l K Hence from the conjecture variant, if \B(k)\ < \B(k+ 1)|, then 
for every u £ B(k), x\u is still in B(k + 1). Therefore it belongs to case (1) we 
reviewed in Section 16.31 and the corresponding column in T\ is a sparse one. If 
\B(k)\ > \B(k+ 1)|, we will have \B(k)\ — \B(k+ 1)| dense columns which come 
from the fact that they belong to case (2) or (3). 

As the coefficients in the the expansion of (1+zH \-z^ d ~ l ^) n are symmetric 

to the central coefficient (or the central two when (d — \)n is odd), the condition 
\B(k)\ > \B(k+ 1)| holds for the first time when k = ko, the index of the central 
term (or of the second one in the central two terms). Then the number of dense 
columns is 

(\B(k )\-\B(k + l)\) + (\B(ko + l)\-\B(k + 2)\) 
+ ... + (\B(k'-\)\-\B(k')\) + \B(k')\ = \B(k Q )\. 

That ends the proof, for such a coefficient |fi(&o)| is exactly the greatest one.D 

The Hilbert series is usually used to analyze the behaviors of Grobner basis 
computation, for example the regularity of the input ideal. As the leading terms 
of polynomials in the Grobner basis and the canonical basis determine each other 
completely, it is also natural to have Proposition 16.41 which links the canonical 
basis and Hilbert series. 
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Remark 6.1 When d = 2, the number of dense columns in T\ is the binomial 
coefficient C*° , where 

n/2 if n is even; 



n+\ 
2 



if n is odd. 



For the case d = 3, such the greatest coefficient is called the central trinomial 
coefficient. 

Corollary 6.5 If the Moreno-Socias conjecture holds, then the percentage of nonzero 
entries in T\ for a generic sequence of degree d is bounded by (niQ + 1)/A where 
mo is the number of dense columns computed from Proposition ^. 4\ 



Proof The number of nonzero entries in the dense columns is bounded by m^D, 
and that in the other columns is smaller than D. □ 

Assuming the correctness of the Moreno-Socias conjecture, we can take a step 
forward from Proposition 16.41 That is, we show case (3) will not occur during the 
construction of 7\ . 

Proposition 6.6 Follow the notations in the Moreno-Socias conjecture. If the con- 
jecture holds, then for any term u G" lt( (P) ), x\ u is either not in lt( (P ) ) or a minimal 
generator oflt((P)). 

Proof Suppose x\U = (mi,. . . ,u n ) G lt((P)) is not a minimal generator. We will 
draw a contradiction by showing u G lt((P)) under such an assumption. 

Without loss of generality, we can assume each u\ 7^ for i = 1 , . . . , n, otherwise 
we can reduce to the n — 1 case by ignoring the z'fh component of u. As x\u is 
not a minimal generator, there exist a&(l <k <n) such that w- k > := («i , . . . , Uk — 
l,...,w„)isinlt((P)). The case when&= 1 is trivial. Otherwise, since deg(ti^) = 
deg(-u) and vS k > < u, by Proposition 16. 3 1 we know u G lt((P}). □ 

Corollary 6.7 If the Moreno-Socias conjecture holds, then the number of dense 
columns in T\ for generic sequences is equal to the cardinality of{g G G\ : x\\ lt(g)}, 
where G\ is the Grobner basis w.r.t. DRL. 

Remark 6.2 By Corollary 16.71 for generic sequences, to construct T\ one only 
needs to find the leading term of which polynomial in Gi is a given term x\u. Thus 
we can conclude that the construction of T\ is free of arithmetic operations. Even 
for real implementations, the cost for constructing T\ is also quite small compared 
with that for the change of ordering (see Section|7]for the timings). Bearing in mind 
that the ideal generated by a generic sequence is in shape position, we know the 
complexity in Theorem l3.2l is indeed the complete complexity for the change of or- 
dering for generic sequences, including construction of T\, the only multiplication 
matrix needed. 
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6.4 Asymptotic analysis 

Next we study the asymptotic behavior of the number of dense columns in T\ for 
a generic sequence of degree d, with n fixed and d increasing to +00. These re- 
sults are mainly derived from a more detailed asymptotic analysis of coefficients 
of the HUbert series of semi-regular systems in [2 3], where standard methods in 
asymptotic analysis, like the saddle-point and coalescent saddle points methods, 
are applied. 

The target of this subsection is to find the dominant term of the greatest coef- 
ficient in the expansion of the Hilbert series H(n,d) in (fTBT ). as d tends to +00 and 
n is fixed. First one writes the with coefficient Id{m) of H(n,d) with the Cauchy 
integration: 

1 j; H(n,d)(z) J 1 J (l-z d ) n 
Id{m) = 2^f ~^^~ dz = 2^f {l-zYz^ dZ - 

With F(z) := / 1-7 yi m+i = e an d g(z) '■= 1, h{m) becomes the form conve- 
nient to the asymptotic analysis 

Idipi) = Y^f s{z)e t{z] dz- 

Suppose the greatest coefficient in H(n,d) comes from the moth term. Since we 
are interested in the asymptotic behavior, we can assume niQ = (d — l)n/2. As a 
special case of [2, Lemma 4.3.1], we have the following result. 

Proposition 6.8 Suppose m$ = (d — l)n/2. Then the dominant term ofld(mo) is 



where f(z) = nlog(^- I T) — (m+ 1) log(z), and ro is the positive real root of f'(z)- 
Furthermore, ro tends to 1 as d increases to +00. 

To prove the fact that the positive real root ro of f'(z) tends to 1, one needs to 
use the equality mo = (d — l)n/2. Other parts of the proof are the same as those in 
HH Section 4.3.2]. 

Next we investigate the value of /"(ro) and F(ro) in the dominant part of 
Id{m Q ) as r tends to 1. Set h(z) := ^^ = 1 +z-\ hz^ -1 . Then 

h"(z)h(z)-h'(z) 2 d+\ 

h(z) z z 

Noting that h{\) = d, h'{\) = d(d - l)/2, and h"{\) = d(d - l)(d-2)/3, we have 

f"(l)=nd 2 / 12 + 0(d). 

With the easily obtained equality F(l) = d n , we have the following asymptotic 
estimation of Id(pio). 
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Corollary 6.9 Let n be fixed. As d tends to +oo, /^(mo) 



6 d"-\ 



This asymptotic estimation of the greatest coefficient in H(n,d) accords with 
the theoretical one. Figure [3] shows the number of dense columns derived from 
both Proposition 16.41 and Corollary 16.91 As can be shown from this figure, the 
asymptotic estimation is good, even when d is small. 



n = 3, d=1..100 



n = 4, d=l.. 100 



ugi/y;) 
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> Theoretical -\- Asymptotic | 



Figure 3: Number of dense columns in 7\ for n = 3,4 and d = 1, . . . , 100 

Corollary 6.10 Let n be fixed. As d tends to +oo, if the Moreno-Socias conjecture 
holds, then the following statements hold: 



1. the percentage of nonzero entries in T\ is ~ </ ^Ljd; 

2. for a generic sequence of degree d, the complexity in Theorem \3.2\ is 0{\ ^D 2+ ~ 



As Corollary I6.10l shows. for a generic sequence, the multiplication matrix T\ 
become sparser as d increases. Furthermore, the complexity of Algorithm is 
smaller in both the exponent and constant compared with FGLM. 

Remark 6.3 Here we only consider the case when n is fixed and d tends to +oo, 
while the asymptotic behaviors of the dual case when d is fixed and n tends to +oo 
have been studied in [3] for the special value d = 2. 



7 Experiments 

The first method for the shape position case, namely Algorithm [2 has been imple- 
mented in C over fields of characteristic and finite fields. A preliminary imple- 
mentation of the BMS-based method for the general case has been done in Magma 
over large finite fields. Benchmarks are used to test the correctness and efficiency 
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of these two methods. All the experiments were made under Scientific Linux OS 
release 5.5 on 8 Intel(R) Xeon(R) CPUs E5420 at 2.50 GHz with 20.55G RAM. 

Table |7] records the timings (in seconds) of our implementations of F5 and Al- 
gorithm |2] applied to benchmarks including theoretical ones like Katsura systems 
(Katsurarc) and randomly generated quadratic polynomial systems of n variables 
(Randoms), and practical ones like MinRank problems from Cryptography lfl6l 
and algebraic cryptanalysis of some curve-based cryptosystem (Edwards). In this 
table, D denotes the degree of the input ideal, and the column "Density" means 
the percentage of nonzero entries in the multiplication matrix T\ . The instances 
marked with f are indeed not in shape position, and the timings for such instances 
only indicate those of computing the univariate polynomial in the LEX Grobner 
basis. The performances of the DRL Grobner basis computation and FGLM in 
Magma (version 2-17-1) and Singular (version 3-1-2), together with the speedup 
factors of our implementation for the change of ordering, are also provided. 

As shown by this table, the current implementation of Algorithm|2]outperforms 
the FGLM implementations in Magma and Singular. Take the Randoml3 instance 
for example, the FGLM implementations in Magma and Singular take 10757.4 
and 19820.2 seconds respectively, while the new implementation only needs 193.5 
seconds. This is around 54 and 101 times faster. Such an efficient implementation 
is now able to manipulate ideals in shape position of degree greater than 40000. It is 
also important to note that with this new algorithm, the time devoted to the change 
of ordering is somehow of the same order of magnitude as the DRL Grobner basis 
computation. 
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Magma 


Singular 


Speedup 
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D 


Density 


Fs(.Q 


New Algorithm 


Ft 


FGLM 


DRL 


FGLM 


Magma 


Singular 


Katsura 11 


2" 


21.53% 


4.9 


3.4 


18.2 


178.6 


632.0 


328.4 


52.7 


96.9 


Katsura 12 


2 12 


21.26% 


31.9 


26.3 


147.9 


1408.1 


5061.8 


2623.5 


53.6 


99.8 


Katsura 13 


2" 


19.86% 


186.3 


189.1 


1037.2 


10895.4 






57.6 




Katsura 14 


2 w 


19.64% 


1838.9 


1487.4 


9599.0 


87131.9 






58.5 
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2 15 


18.52% 


11456.3 


12109.2 














Random 1 1 


2" 


21.53% 
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3.4 


18.1 


169.3 
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328.6 


49.2 


95.5 
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2 12 


21.26% 


26.6 
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134.9 


1335.8 


4867.4 


2581.1 


49.6 


95.8 
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2" 


19.98% 


146.8 


193.5 


949.6 


10757.4 


36727.0 


19820.2 


55.6 


102.4 
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2"> 


19.64% 


1000.7 


1489.5 


7832.4 


84374.6 






56.6 
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6882.5 


10914.02 
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1.1 


0.5 


6.3 


22.7 
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38.1 


43.6 
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4116 


22.95% 


28.4 


28.5 


208.1 


1360.4 


4985.8 


2490.3 


47.7 


87.4 


MinR(9,8,5) 


14112 


19.04% 


543.6 


1032.8 
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41580 


16.91% 


9048.2 


22171.3 
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4096 


7.54% 


4.0 
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5.8 


418.3 


72.4 


1823.6 


46.7 


203.7 


Edwards \ 


4096 


3.41% 


0.1 


2.4 


0.2 


176.7 


1.0 


839.9 


72.7 


345.6 


Cyclic 10 f 


34940 


1.00% 




3586.9 


>16hrs 


and >16 Gig 











Table 2: Timings of the method for the shape position case from DRL to LEX 

Table [3] illustrates the performances of Algorithm 5] for the general case. As 
currently this method is only implemented preliminarily in Magma, only the num- 
ber of field multiplications and other critical parameters are recorded, instead of 
the timings. 

Benchmarks derived from Cyclic 5 and 6 instances are used. Instances with 
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ideals in shape position (marked with J) are also tested to demonstrate the general- 
ity of this method. Besides n and D denoting the number of variables and degree of 
the input ideal, the columns "Mat Density" and "Poly Density" denote the maximal 
percentage of nonzero entries in the matrices 7\ , . . . , T n and the density of result- 
ing Grobner bases respectively. The following 4 columns record the numbers of 
passes in the main loop of Algorithm |4j matrix multiplications, reductions and field 
multiplications. 

As one can see from this table, the numbers of passes accord with the bound 
derived in theorem I4TT1 and the number of operations is less than the original FGLM 
algorithm for Cyclic-like benchmarks. However, for instances of ideals in shape 
position, this method works but the complexity is not satisfactory. This is mainly 
because the resulting Grobner bases in these cases are no longer sparse, and thus 
the reduction step becomes complex. Fortunately, in the top-level algorithm |5J it is 
not common to handle such ideals in shape position with this method. 



Name 


n 


D 


Mat Density Poly Density 
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#Mat. 


#Red. 


#Multi. 
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nD 2M 
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65 


8.73% 


19.7% 
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nD 2.m 
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4 


36 


60.65% 
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nD 2m 


Dl± 


12 


48 


34.2% 


51.02% 


624 


780 


517 


mD 2.874 


Dessin2-6 \ 


6 


42 


46.94% 


100% 


294 


336 


205 


nD 2968 



Table 3: Performances of Algorithm |4]for the general case from DRL to LEX 
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